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ABSTRACT 

We investigate the power spectra of outflow-driven turbulence through high-resolution 
three-dimensional isothermal numerical simulations where the turbulence is driven lo¬ 
cally in real-space by a simple spherical outflow model. The resulting turbulent flow 
saturates at an average Mach number of ^ 2.5 and is analysed through density and 
velocity power spectra, including an investigation of the evolution of the solenoidal 
and compressional components. We obtain a shallow density power spectrum with a 
slope of ^ —1.2 attributed to the presence of a network of localised dense filamentary 
structures formed by strong shock interactions. The total velocity power spectrum 
slope is found to be —2.0, representative of Burgers shock dominated turbulence 
model. The density weighted velocity power spectrum slope is measured as ^ —1.6, 
slightly less than the expected Kolmogorov scaling value (slope of -5/3) found in pre¬ 
vious works. The discrepancy may be caused by the nature of our real space driving 
model and we suggest there is no universal scaling law for supersonic compressible 
turbulence. We find that on average, solenoidal modes slightly dominate in our turbu¬ 
lence model as the interaction between strong curved compressible shocks generates 
solenoidal modes, and compressible modes decay faster. 
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1 INTRODUCTION 

Observations, theory, and simulations have shown that the 
physics of star forming regions are dominated by supersonic 
turbulence. The turbulence supports against gravitational 
collapse of the molecular cloud and leads to a hierarchi¬ 
cal structure of diffuse regions, filaments, and clumps, as 
kinetic energy is distributed fro m kilo-parsec scales t o sub- 
Astronomical Unit (AU) scales jBurkhart et al.11^131 ). Sim¬ 
ulations have shown supersonic turbulence to decay rapidly, 
in less than a crossing time, and so the turbulence must be 
continuously replenished in order to account for the rela¬ 
tively low star fo rmation rates (SFR) that are observed in 
molecular clouds (iMac Low fc Klessenll20(3 ). 

Protostellar jets and outflows, being synonymous with 
star formation, are one possible candidate as a mechanism 
to drive and sustain supersonic turbulence within molec¬ 
ular clouds on small scales. They are an essential compo¬ 
nent of the early stages of star formation. As the collaps¬ 
ing protostellar system expels an excess of angular momen- 
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turn in the form of bipolar protostellar jets along the ro¬ 
tation axis of the system, the jets entrain and accelerate 
ambient cloud material to form the protostellar outflows. 
The outflows emit profusely in molecular lines revealing to 
observers information about the local amb ient cloud envi¬ 
ronm ent and the embedded driving sources llReiter fc SmithI 
I2OI3I) . Although observations suggest outflows may have 
little effect on the largest scale turbulence in molecular 
clouds (e.g. iMa c Low fc Klesseiil l2004l: iPadoan et al.l l2009l : 


IZrake fc MacFadvenI 201? Reiter fc SmithI 20131) . the com¬ 

pressive modes they generate may sti ll be important at trig - 
gering star formation on small scales llSchneider et al.ll201ll ). 


In a study bv iFederrath fc KlessenI ||2012|) . the authors 
found that both Mach number and turbulent forcing are the 
main parameters in controlling the SFR. They establish that 
turbulent compressive type forcing can increase the SFR 
up to an order of magnitude compared to solenoidal type 
forcing, and a Mach number change from 5 to 50 can increase 
the SFR by a factor of up to five. The effects of magnetic field 
were revealed to play a more minor role by only reducing 
the SFR by a factor of two. Although playing a minor role, 
recent studies of models and simulations have shown that 
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the presence of magnetic fields are still required in order 
to redu ce the SFR rate to the levels observed in molecular 
clouds (jPadoan et alJlioi^ l. 

NGC 1333 is an extremely rich star forming region. 
Part of the Perseus molecular cloud complex, it is classed 
as a reflection nebula, covering an area 7° x 3° in the sky, 
and esti mated through VL BI methods to be 235±18 parsec 
distant llHirota et al.11200^ . Its vicinity and diversity allow 
observers to study all phases of the star for mation process 
inclu ding, for example, protostellar collapse l|Belloche et al.l 
l2006l ) , and the accretio n disks and jets of a protobinary sys¬ 
tem dChoi et al.ll201ll l. Mapping observations through CO 
(J = 3 — 2) and (J = 1 — 0) have shown the structure of 
NGC 1333 to be dominated by many molecular outflows that 
have broken the cloud into low-den sity cavities separated by 
dense shells (iKnee fc Sandellll 2000l l. A recent maser survey 
of NGC 1333 by Lvo et al. (2014) found several masers with 
no associated protostar or outflow detected suggesting they 
may indicate a very early evolutionary stage of star forma¬ 
tion before the onset of the traditional outflow phase. Can 
the energy deposited by the outfl ows in NGC 133 3 maintain 
the observed level of turbulence? iMatzneJ (l2007l 'l proposed 
that the observed level of turbulence in NGC 1333 could 
be d riven internally by outfl ows if certain conditions were 
met. iGutermuth et alj ll2008l 'l found evidence that Class II 
sources have dispersed surrounding molecular gas on short 
timescales suggesting outfl ows can ultim a tely d estroy the 
NGC 1333 cloud. However, IPadoan et al.l ll2009h measured 
the turbulent energy spectrum of NGC 1333 ^^CO maps 
and found slopes of ~ —1.85 with no flattening near the 
scale expected for outflow driving, suggesting turbulence in 
NGC 1333 may not be driven by ou tflows. A JCMT survey 
of NGC 1333 bv ICurtis et al.l (|201fll f of C^®0 data estimates 
that there is more energy introduced by outflows than is ob¬ 
served in the bulk cloud, implying that it is not easy to feed 
high-energy outflow gas into turbulent motions. 

In order to understand momentum transfer and obser¬ 
vational power spectra in more detail, we can study it most 
readily and precisely through numerical simulations of tur¬ 
bulence. In such simulations the tur bulence is either driven 
globally through Fourier space (e.g. Federrath et al. I I2OO8I: 
iMicic et al.ll20l3 : iFederrath fc KIessenll2013h. or locally in 
real-s p ace via protostellar outflo w s (e.g. Nakamma&Ml 
2OO7I: Carroll et al.l l 200 ^ . 120101 : iMoraghan et al. 20131 : 


Federrath et al.ll2014h . In this paper we drive the turbu¬ 


lence in_j;eal^ac£with outflows using a parameter scheme 
from iMatzneil ll2007ll based on the observational properties 
of NGC 1333. 

Turbulence consists of eddies that transfer energy across 
a range of scales from the injection scale to the dissipative 
range whilst conserving momentum. The degree of turbu¬ 
lence is dependent on the Reynolds number which describes 
the r atio of the inertial to viscous forces dBurkhart et al.l 
l2013l ~). T he turbulence pheno menon was first statistically de¬ 
rived by iKolmogorovI (Il94ll i through the assumption that 
the velocity field is a stochastic distribution with a constant 
energy transfer rate at all scales. The result is the famous 
Kolmogorov energy spectrum for incompressible turbulence 


E{k) oc 


( 1 ) 


where the energy power spectrum, E{k), is proportional 
to the energy input rate per mass e, and the wavenum¬ 


ber, k (also known as the inverse length scale; k = 
27r/A), always possesses a nega tive power-law index of -5/3 
(iFalceta-Goncalves et al]l2014l i. 

However, within the astrophysical regime - such as 
molecular clouds - turbulence is highly supersonic and com¬ 
pressible and the measured energy power spectra have been 
shown to clearly deviate from the incompressible -5/3 value. 
It was considered that the physics may be more appropri¬ 
ately represented by the Burgers turbulence model. The 
Burgers model is related to a fully analytical description of 
ID shocks in a 1948 paper by Burgers as an exact solution to 
the Navier Stokes equations. The fundamental difference of 
the Burgers turbulence model to Kolmogorov turbulence is 
that energy is transferred from large scales to s mall scales via 
shock s rather than a cascade through eddies llVestuto et al.l 
I2OO3II . It envisages many shocks possessing step functions, 
and their Fourier transform is proportional to k~^ and the 
resulting velocity power spectrum is 

P^{k)(xk-^. ( 2 ) 

In recent years, much work has focused on finding a 
universal scaling law to represent supersonic compressible 
turbulence in the same way as the Kolmogorov law closely 
defines incompressible turbulence. A power spectrum based 
on the density is not suitable as simulations have shown 
a dependence of d ensity with Mach number. For example, 
I Kim RvrJ ([200^ performed isothermal Fourier driven tur¬ 
bulence simulations driven by solenoidal modes for Mach 
numbers in the range of 1 to 12. They found the density 
power spectral slope flattened with increasing Mach num¬ 
ber due to the development of sheet-like and filamentary 
density structures ranging from a slope of -1.73 for A4rms 
= 1.2, to a slope of -0.52 for AIrms = 12. This shows that 
power index for density is a function of Mach number where 
at transonic velocities the density power spectrum may fol¬ 
low Kolmogorov turbulence. Velocity power spectra may 
show a similar trend. Although modern supersonic com¬ 
pressible turbulence simulations typically possess velocity 
power-law exponents of -2 as the turbulence would be shock 
dominated and r epresentative of Burgers turbulence such as 
iFederratiil ll2013li . early simulations of subsonic and tran¬ 
sonic fl£ws_werefitted with spectral indices similar to -5/3 
(e.g. iPorter et al.lll99'3j . 

A highly supersonic velocity power spectral index is 
expected to be close to -2. However, a turbulent veloc¬ 
ity domain is a combination of compressible (dilational) 
mod es, and solenoidal (lon gitudinal or transverse) modes 
(e.g. IFederrath et al.l l2010l h Different driving mechanisms 
consist of different strengths and ratios of the constituent 
solenoidal and compressible components. It has been found 
to have a strong effect on the resulting density distribution 
with compressible driving leading to a broader density distri¬ 
bution i n comparison to one g enerated by purely solenoidal 
driving (IFederrath et al.ll2008lf . This is because compressive 
forcing generates shocks leading to stronger compression and 
raref action in the computational domain (IFederrath et al.l 
2OIOI). It also leads to d ifferent spectral slopes. For example, 
Federrath et al.l (l2010l j obtained a velocity power spectral 
slope of -1.86 for pure solenoidal driving and -1.94 f o r pure 
compressional driving. Further study in iFederraHil (l2013fl 
showed it is difficult to directly compare the power spectrum 
of solenoidal driving to compressible driving due to differ- 
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ent inertial ranges and different dependencies on numerical 
resolution. Due to the spectral differences of the underlying 
modes, a pure velocity spectrum could not be considered as 
a universal scaling relation. 

An alternative relation is being studied as a poten¬ 
tial universal scaling relation; the density weighted velocity 


power spectrum It was first tested bv lKritsuk et al.l 

(l2007l i who found it to scale as in-line with the 

Kolmogorov law. But later simulations showed that there 
could still be a dependence o n the driving mechanism . Fully 
compressible simulations by Federrath_elLal] (|201(]|) found 
P{p^^^v) oc Later. ICaltier &: Baneri^ (l201ll ') theo¬ 

retically predicted that a strong compressible component 
would produce a P{p^^^v) scaling relation with a slope of 
— 19/9, which was in good agreement to the simulations. 
The reason solenoidal driving would provide a different re- 
lation is explain ed by the previously mentioned analysis of 
iFederra^ (120131 ') regarding the inertial-range scaling of the 
compressible components. In this paper we will examine the 
scaling using our real-space outflow driven turbulence 
model. 

Unlike the majority of other numerical based turbu¬ 
lence studies which drive the turb ulence on large global 
scales through Fourier-space (e.g iFederrath et al.l l2010l : 


iKritsuk et alll2013l : iFederrathl l2013l : Yoo fc Ch3 2014l ~) , we 
drive the turbulence exclusively in real-space on small local 
scales with spherical outflows. Fourier-space driven turbu¬ 
lence has close ties to power spectra as the power spectra 
are themselves derived through Fourier space. Our real-space 
model is not only more realistic physically, but also com¬ 
pletely impartial to the power spectra analysis, and thus pro¬ 
vides a useful comparison to results of Fourier-space driven 
turbulence. Our aim is to quantify our real-space outflow 
model and see how it compares to Fourier driven turbulence 
models in terms of power spectra. 

In reality, turbu lence may be driven at more than one 
scale. Simulations bv ICarroll et al.l ll2010lj included an exter¬ 
nal isotropically driven turbulence in addition to their out¬ 
flow driven turbulence model, and obtained steeper velocity 
power spectra with a detectable kn ee in the power spe ctra at 
the outflow driving scale. Recently. I Yoo &: Chd ll2014h stud¬ 
ied Fourier driven turbulence driven at two separate scales 
simultaneously, effectively accounting for different driving 
mechanisms such as supernovae or galactic rotation shear 
on larger 100 pc scales, and stellar winds and protostellar 
outflows on smaller parsec scales. They found the effects of 
the smaller scale driving may not be easily observed through 
the kinetic energy spectrum unless it is of comparable en¬ 
ergy to the larger scale driving mechanism. However, due to 
our localised driving, our method can only be used to as¬ 
sume small scale local turbulent driving near a star-forming 
clump within a molecular cloud. 

Our model produces an average density power spec¬ 
trum slope of —1.2, a lower value than is usually found 
in other works. On the contrary, we measure the average 
velocity power spectrum slope to be -2.0, a value agreeable 
with that of Burgers shock dominated turbulence. When we 
plot the density weighted velocity power spectrum 
we obtain a slope of -1.6. This value is slightly shallowe r 
than the pure solenoidal driving case of IFederrathl ll2013lj . 
an d significantly shall o wer t han the value derived 

by lOaltier k. Baneried ([iolj), and obtained through simu¬ 


lations by IFederrathl (l2013fl for purely compressive driving. 
It suggests a fundamental difference between Fourier based 
driving and real-space driving. However, our measured slope 
of -1.6 may be comparable to the Kolmogorov scaling value 
of -5/3 (~ 1.67). The small difference in the measurement 
may not be significant given the uncertainties from time 
variations and systematics introduced by our relatively low 
resolution and resulting short inertial range for fitting the 
slope. 

In terms of compressional and solenoidal components, 
we find strong compressional modes located at active out¬ 
flow bowshocks and instances of outflow-outflow collision. 
Solenoidal modes are also generated at these locations, 
mostly through curved shocks, before dispersing through¬ 
out the computational domain. Overall we find compres¬ 
sional modes have a slightly shallower power-law slope, sug¬ 
gesting these modes are more efficient at transferring en¬ 
ergy from large scales to small scales. This is supported 
thro ugh the ratio of the so lenoidal to compressional modes 
(^See IFederrath et al.ll20(I^ . which shows solenoidal modes 
are created through our outflow-outflow shock interactions 
and ultimately slightly dominate on average throughout the 
computational domain in an attempt by the system to reach 
an equilibrium based on the degrees of freedom of the differ¬ 
ent modes. However, at instantaneous times, both solenoidal 
and compressional components can be highly variable and 
interchangeable during the actively driven simulation. 

The layout of our paper is as follows; We describe our 
computational code and model setup in § [51 numerical re¬ 
sults and power spectra are presented and discussed in § (3) 
and we conclude in § |T] 


2 NUMERICAL METHODS 

2.1 Numerical code and basic equations 

The simulations presented here were performed using an 
isothermal multi-dimensional fixed-grid code based on the 
total variation diminishing scheme, known as the tvd code. 
It is an explicit finite-difference code and uses a second-order 
accurate Roe-typ e up-winded Rie mann solver to calculate 
the inter-cell flux(p<irn et al.lll999ll . 

iKitsionas et al.l ll2009ll compared the performance of the 
TVD code against three contemporary grid-based codes - 
ENZO, FLASH, and ZEUS - through simulations of decaying, 
isothermal, supersonic turbulence where each code used the 
same initial conditions. It was found that the tvd code pro¬ 
duced comparable results to the others, and more efficiently 
with shorter run-times. This efficiency makes the tvd code 
particularly suitable for our outflow driven turbulence sim¬ 
ulations. As there are a wide range of velocities on the grid 
at any given time, from the supersonic motion of the active 
outflows to the subsonic motions of decaying turbulence in 
remoter parts of the grid, the computational time-step is 
constrained by the fastest motions and the simulation must 
also be run for many dynamical timescales. Therefore use of 
a computationally efficient code is preferable. 

We use an isothermal version of the tvd code. One may 
expect dynamical processes in molecular clouds to include 
molecular cooling a nd chemistry. However, a study of molec¬ 
ular turbulence by IPavlovski et al.l (l2006l l found that the 
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isothermal approximation can be sufficient as the cooling 
timescale is shorter than the dynamical timescale. 

Previous works concerning numerical t urbulence simu¬ 
lations have included magnetic fields (e. g. iLiJ^^jakamural 
200^ Wang et"^ 2010l : ICho fc Kimll201ll : Yoo fc Ch^ 20 14 


Federrath et al.ll2014l '). In this paper, we however only show 
spectral characteristics of a hydrodynamic turbulent flow 
driven by outflows in real-space. In a later study, we will 
include magnetic fields. 

The isothermal tvd code solves the following two equa¬ 
tions of hydrodynamics presented in Euler conservation 
form. The conservation of mass, 


dtp -I- V • (pu) = 0 (3) 

and the conservation of momentum, 

dtv + V ■ Vu = Cs + F (4) 

P 

where p is density, t is time, Ca is the isothermal sound speed, 
V is the flow velocity vector, and F is the outflow forcing 
term. 

Our primary simulation is performed on a fixed uni¬ 
form grid of 1024^ cells with periodic boundary conditions. 
This effectively emulates a small sub-set of a larger molec¬ 
ular cloud. We repeated the simulation at lower resolutions 
of 512® and 256® cells for a convergence test as resolution is 
an important factor when analysing data via power spec- 
tra in order to resol ve an inertial range sufhciently (e.g. 
iFederrath et al.lf2010ll . 


zero and the total momentum on the computational domain 
remains conserved. 


2.3 Dimensionless Parameter scheme 


We implement the same dimens io nless parameter 
sc heme first intro d uced by iMatzneil | 2003 ) and used 
in iMoraghan et al. I ll2013ll to provide a phy sical measure 
to our mode l. This scheme was also used by I Carroll et al.l 
ll20n9l. 1201(1) and ties our results to physical quantities. 
It is in fact based on the physical parameters measured 
from NGC 1333. The dimensionless quantities of mass, mo, 
length, lo, and time, to, are defined as follows; 


4„3 3 

_ p7P7 Pt , , _ 

TTIQ — o 5 tQ 11? &Ild. 60 3 4 5 


( 6 ) 


with parameters of density, outflow momentum, and outflow 
rate per unit volume set as p = 2.51 x 10“^® gcm“®, P = 
3.98 X 10®® gcm“®s“®, and S — 6.31 x 10~®® cm“®s“®, 
respectively. The resulting physical unit quantities we use 
are thus; mo = 18.7 Mq, Iq — 0.37 pc, to = 0.34 Myr, and 
Vo = lo/to ~ 1.064 kms“®. W ith these values the loca l sound 
speed is 0.587 kms“®. As in iMoraghan et ^ ll2013li we set 
each side of the computational box as L = 8Zo (3.96 pc), 
but this time we choose an outflow radius of P = 0.35 lo (a 
lower limit necessary for numerical stability of our code at 
1024®), and found it sufficient to run the simulation for time 
2.0 to (-0.68 Myr). 


2.2 Forcing scheme 

Our turbulence forcing model proceeds as follows; A 3 di¬ 
mensional (3D) computational domain is initialised with a 
uniform gaseous medium. Different points within the domain 
are randomly selected. At each randomly selected point, the 
mean density, p, of a volume defined by the outflow radius, 
R, is determined. If the mean density of this volume is larger 
than the global mean, po, and no point within the volume is 
below a lower-threshold value O.OOlpo, then it is chosen to 
be the location of an outflow. 

The momentum, Pj, which an outflow should introduce 
is defined as follows, 

III 

<R 

where p(f) is the density within the outflow volume depend¬ 
ing on the distance vector r from the outflow center, and 
Vj (P) is the outflow velocity at its outermost boundary, P. 
The exponential q parameter, set to 1.0, attempts to mimic 
the observed ‘Hubbl e-law’ type veloc ity observed in proto- 
stellar outflows fe.g. lArce et al.l[2007h . 

As the total scalar momentum, P, which each outflow 
should introduce is known, Vj (P) can be expressed in terms 
of P and the integral part of Eq.[5] However, by setting a 
new velocity field, the total momentum of the outflow may 
not be conserved due to a non-uniform density distribution 
existing within the chosen outflow radius, especially dur¬ 
ing later times of the simulation. Therefore we calculate the 
added momentum and then subtract it within the radius 
P. This ensures the vector sum of the added momentum is 


2.4 Power Spectra 

In the next section we will analyse our turbulent data in 
terms of time-averaged compensated power spectra. The 
power spectra are time averaged over a total of 50 output 
dumps between time 1.0 - 2.0 to, a range where the turbu¬ 
lence has fully saturated whilst being actively driven. The 
compensated po wer spectra method is often used in turbu¬ 
lence studies ( e.g.lKitsionas et al.| l2009l : IZrake fc MacFadvenI 
I 2 OI 2 I : IFederrath fc Klessenll2013ll . as it aids in visualising the 

correct spectral index and inertial range. The power spec¬ 
trum of a simulated quantity, P{k), is compensated by a 
factor, , where the exponent S is chosen to negate the 
power-law slope, or in other words, the power-law exponent 
is the negative of 5. We only use the highest resolution sim¬ 
ulation to determine the power-law slope. Due to the nature 
of our real-space driving occurring on small scales, an ap¬ 
propriate inertial range is quite limited in comparison to 
Fourier based driving. 

The three-dimensional velocity field can hold additional 
information as to the nature of the velocity motions. In real 
space a velocity field can be decomposed into three direc¬ 
tions; X, y, and 2 ;. Correspondingly, in Fourier space, it can 
be decomposed into two; solenoidal (also known as ‘diver¬ 
gence free’), and compressible (or ‘curl free’). The mathe¬ 
matical principle is based on the Helmholtz decomposition 
theory which states that a differentiable vector held with 
divergence (V • u) , and curl [V x u], can be split into its 
corresponding compressional and solenoidal components 

v{r) = Vs(r) + Vc{r) (7) 

where V x Uc = 0 and V • Us = 0. The resulting solenoidal 
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Time [to] 

Figure 1. Root-mean-square Mach number evolution with time 
during our simulation at all resolutions; 1024® (black), 512® 
(blue), and 256® (red). The turbulence saturates after a time of 
0.5 to- The vertical lines at 1.0 to and 2.0 to indicate the time inter¬ 
vals between which we consider our time-averaged power spectra. 
The average RMS Mach number between time 1.0-2.0 to is mea¬ 
sured as 2.456 for the 1024® simulation, 2.596 for the 512®, and 
2.460 for the 256® 

and compressional power spectra are defined as, 



(8) 


(9) 


where Vs{k) and Vc{k) are the Fourier transforms of the 
solenoidal and compressional components, respectively. Vor- 
ticity in the turbulent flow leads to solenoidal components, 
and converging flows generate compressible components. 
Studying these components and their ratios is a method 
now regularly used to characterise and understand turbu¬ 
lence in numerical simulations (e.g. iFederrath et al.l l2008l . 
I 2 OI 0 I : iMicic et al.ll2012l : lFederrathll2013l ). 

3 RESULTS 

3.1 RMS evolution and 2D visualisations 

Our simulation quickly generates and maintains supersonic 
turbulence at a saturated level. This can be seen in the plot 
of root-mean-square (RMS) Mach number (AIrms) evolu¬ 
tion with time in Figure[T] For comparison, the AIrms evo¬ 
lution of the 512^ and 256® simulations are plotted as blue 
and red lines, respectively. In all cases supersonic turbulent 
motions become fully saturated after a time of 0.5 to. We run 
the simulations to a time of 2.0to, at which point we have 
obtained 50 output dumps of fully saturated turbulence be¬ 
tween a time of 1.0-2.0 to for use in the time averaged power 
spectra analysis to follow. The high-velocity spikes are a con¬ 
sequence of our real-space driving model. When an outflow 
is added to the grid, a small region may consist of very low- 
density cells. Through the conservation of momentum, these 
cells would have a much higher velocity than average leading 
to a high velocity spike. However, the high velocity cells are 
highly localised and do not effect the global turbulence. 

To aid visualisation of the data, a slice through the 3D 
data cube representing the log of density at time 2.0 to is 
shown in Figure[2] Here we see the turbulent filamentary 



0 2 4 6 8 

^ [Iq] 


Figure 2. Slice through the 3D volume depicting log of density at 
time 2.0 to for the 1024® simulation. The outflows sweep the am¬ 
bient gas into thin shells which interact with other shells creating 
the filamentary structure interspersed with lower density voids. 
Arrows represent the velocity field vectors for high velocities. 

structure as well as newly expanding low-density outflow 
cavities. The expansion of the outflow cavities are high¬ 
lighted by the arrows representing the velocity field on the 
displayed 2D plane for high velocities. These velocity vectors 
show the velocity profile of the outflows range from zero ve¬ 
locity at the interior, to maximum velocity at the boundary. 
We see that high velocities are not sustained far from the 
outflow sources suggesting most of the momentum is lost 
through the bow-shock and converted into slowly expand¬ 
ing shells. This phenomenon is more clearly seen in Figure|3J 
which plots the divergence (V ■ v) of the velocity field. The 
divergence field highlights expansion flows as positive values 
(light colour) and areas of compression as large negative val¬ 
ues (dark colour). We see the strongest shocks occur around 
the boundaries of the new outflow cavities which are im¬ 
mersed in a filamentary structure of older weaker shocks. 
The weake r shock structures are rem iniscent of the ‘fossil 
cavities’ of ICunningham et al.l ll2006l) . In a complimentary 
fashion, we can plot the curl (V x v) of the velocity field, 
where the curl represents rotational motions. Figure|4] shows 
the 2 -component of the curl V x u revealing that whilst most 
of the domain is saturated by small absolute values of ro¬ 
tational motion, interacting regions possess extreme values. 


3.2 Compensated power spectra 

Now we analyse the full data cubes in terms of power spec¬ 
tra. We plot the compensated power spectrum of density in 
Figure[3 and velocity in Figurej^ 

For our density power spectrum, Pp{k), a slope of -1.2 
was found to be the best fit to the data of the 1024® sim¬ 
ulation between an inertial range of fc = 20 — 40 as seen in 
FigureO The slope of the power spectrum flattens below this 
scale as usually seen in power spectra of outflow driven tur- 































6 


A. Moraghan, Jongsoo Kim & Suk-Jin Yoon 




^ [Iq] 

Figure 3. Slice through the 3D volume depicting real-space diver¬ 
gence (V-u) data at time 2.0 £o for the 1024^ simulation. Expand¬ 
ing flows have positive divergence (light colour), whereas converg¬ 
ing flows have negative divergence (dark colour). This highlights 
the fact that the strongest shocks occur at the bow-shocks of 
expanding outflow cavities. 


/ ' 
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Figure 4. Slice through the 3D volume depicting the z- 
component of the curl, V x u, at time 2.0 to for the 1024^ sim¬ 
ulation. This plot highlights the intricate structure of rotational 
motions which have been generated within the entire volume. 

bulence sim ulations unless an additional larger scale driver is 
present fe.g. lYoo fc Chdl2014l) . The vertical dashed line rep¬ 
resents the outflow driving wave-number, koutfiow/ko. This 
is the scale on which the outflows introduce their energy. 
It is deflned as the outflow size expressed in wavenumber 
koutfiow = 2tt/R, divided by the size of the computational 
box expressed in wavenumber, ko = 27r/L. We also plot the 
results of the lower resolution 512® (blue) and 256® (red) 
simulations. With decreasing resolution there is a decrease 
in the number of high-frequency wave-numbers, limiting the 


Figure 5. Compensated power spectrum of the time averaged 
density data for the 1024® simulation (black). The normal power- 
spectrum, Pp(fc), is compensated by the factor where S = 1.2, 
represents the ideal power-law slope taken at an inertial range 
between fc = 20 — 40. The corresponding compensated density 
power spectra for the 512® (blue) and 256® (red) simulations are 
also shown at the same value of . A suitable inertial range 
is quite limited as our real-space driving method occurs on small 
scales only. The vertical dotted line represents the outflow driving 
scale. 


suitable inertial range and showing the importance of high- 
reso lution simulatio ns for accurate power spectral analysis. 

iMatzneJ ll2007l) predicts the spectral slope of density 
would increase in regions where outflow momentum is cap¬ 
tured more effectively. The many low-density cavities cre¬ 
ated by the strong shocks of the outflows on our grid may 
translate to a lower than expected power spectrum. A large 
volume of low-density gas in real-space translates as an ab¬ 
sence of lower frequency wave-numbers in Fourier space and 

hence a shallowing of the spectral slope. _ 

Our result is intermediate to that of iKim fc Rvul (l2005l) 
whose Fourier driven simulations obtained density spectral 
slopes of -1.73 for AIrms = 1-2 turbulence, to a flatter -0.52 
for AIrms = 12.0. The authors accounted for this flattening 
as a consequence of the formation of dense filaments and 
sheets as Mach number increases. If we assume a smooth re¬ 
lationship bet ween spectral slop e and Mach number in the 
simulations of iKim fc Rvul (l2005l ) , we can perform interpola¬ 
tion to determine the spectral slope their model should pro¬ 
vide if their turbulence was driven at the same Mach number 
as our model. After interpolation using a cubic polynomial 
fitted to their data points we determine that a A^rms = 
2.456 flow would lead to a spectral slope of -1.3. This is 
slightly steeper than our value of -1.2 for the same Mach 
number. Several factors may account for this difference (e.g. 
driving mechanism and chosen inertial range). 

The compensated power spectrum of velocity is pre¬ 
sented in Figure|6] Here we find the power spectrum slope 
for the 1024® simulation is best compensated by a 5 ex¬ 
ponent of 2.0 for an inertial range of fc = 20 — 40. This 
time there is greater deviation between the spectral index at 
different resolutions. The lower-resolutions appear to show 
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Figure 6. Compensated power spectra of the averaged time ve¬ 
locity data of the 1024® (black), 512® (blue), and 256® (red) sim¬ 
ulations. In a similar fashion to Figure [S] the power-spectrum 
of velocity, P„(fc), is compensated by the factor P where 3 rep¬ 
resents the ideal power-law slope between an inertial range of 
fc = 20 — 40 for the 1024® data, which we find to be 5 = 2.0. 
At lower-resolutions, the velocity power spectra display poorer 
convergence compared to the density power spectra. The vertical 
dotted line represents the outflow driving scale. 


more power between the inertial ranges of fc = 10 — 30 . 
This is also noticeable in the Figure 6 of iFederrathI (l2013l ) 
who analysed power spectra of several resolutions up to an 
extremely high-resolution of 4096^. However, our compen¬ 
sated spectral slope is comparable within the chosen, albeit 
limited, inertial range of fc = 20 — 40 at our highest reso¬ 
lution as the 1024® resolution provides no suitable inertial 
range below fc = 20. Overall the value of -2.0 for the slope 
suggests the velocity spectrum represents that of Burgers 
type turbulence and hence is dominated by strong shocks. 
This value of spectral index of velocity is frequently found 
in many previous studies at higher Mach num bers. It was 
obtained, for example, by Kritsuk et ^ (l2007l ) with Mach 
number 6 driving, and IFederrathI ( 20131) with Mach number 
17 driving. In comparison, our RMS Mach number is rel¬ 
atively low at 2.456. This suggests that the total velocity 
power spectrum is not dependent on Mach number once a 
turbulent flow is supersonic. 



Figure 7. Average density-weighted velocity power spectrum 
(p®/®t!) of the 1024® simulation. We fit a slope of -1.6 between the 
inertial range of fc = 20 — 40 for the 1024® simulation. The result 
of the 512® (blue) and 256® (red) simulations are also shown. The 
vertical dotted line represents the outflow driving scale. 


the driving scale possibly through compressible shocks cross¬ 
ing multiple scales, unlike shocks from solenoidal driving. 
Therefore it is highly dependent on numerical resolution. 

When we plot the average density weighted velocity 
power spectrum of the 1024® simulation in Figure!?] we 
obtain a spectral slope of -1.6 in the previously chosen iner¬ 
tial range of fc = 20 — 40. We may have expected our result to 
be between the fc"^'^"* pure soleno i dal an d fc“^'^® pure com- 
pressional extremes of iFederraHtl ll2013l ) as our real-space 
driving is a combination of both solenoidal and compres- 
sional modes. But our result is actually shallower than their 
pure solenoidal limit. In addition, our result is slightly shal- 
lo wer than the -5 / 3 ('^ 1.67) scaling as was initially expected 
bv iKritsuk et al.l (l2007l) . although the difference may not be 
significant given the short inertial fit range and uncertainties 
arising from possible time variations and resolution effects of 
our model. As our localised real-space driving model mostly 
agrees with Fourier based driving through the velocity field, 
it is likely the deviation in P(p®'^®ii) could be due to our den¬ 
sity field where either our real-space localised driving may 
create a different density field morphology than one based on 
global Fourier driving, or the densi ty field is very dep endent 
on Mach number as was found in dKim fc Rvull20n^ ). 


3.3 Density weighted velocity power spectra 

Density weighted velocity power spectra (p^/®n) have been 
proposed as a way for the Kolmogorov law to hold for com¬ 
pressible turbulence due to the mixed quantity statistics in¬ 
volved. A detailed study of densit y weighted velocit y power 
spectra was recently performed bv IFederrathI ll2013l) and re¬ 
vealed it may not be a universal scaling relation as first 
thought. Although it is expected that P(p®/®u) oc fc-®/® for 
low resolutions, the authors found a density weighted veloc¬ 
ity power spectrum of for pure solenoidal driving, and 

a steeper for pure compressional driving. The reason 

suggested by the author for this difference is that compress¬ 
ible turbulence has a wider inertial range and depends on 


3.4 Solenoidal and compressional velocity 
components 

We can also compare our turbulence model to that of Fourier 
driven turbulence in more detail by separating the velocity 
field into the corresponding solenoidal and compressional 
components. The total velocity power spectrum is presented 
in Figure|8]in addition to its constituent solenoidal (red) and 
compressional (blue) components. Once again, we see a slope 
of ~-2.0 fits the total velocity power-spectrum between the 
inertial range of fc = 20 — 40. Our compressible modes slope 
of -2.05 is slightly steeper than the s olenoidal modes slop e 
of -1.98. This trend was also found in iKritsuk et al.l ll2007h . 
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Figure 8. Averaged total velocity power spectrum (black line) 
and its solenoidal (blue) and compressional (red) components for 
the 1024^ simulation. Vertical dotted line represents the outflow 
driving scale. We average velocity data of 50 output dumps be¬ 
tween time 1.0-2.0 to- Power-spectra slopes measured over an in¬ 
ertial range of fc = 20 — 40. We see a slight increase of power 
in the low wavenumbers above the outflow driving scale due to 
the physical expansion of the outflows. They effectively transfer 
power from their initial driving scale to a larger scale. This effect 
is not detectable through the compensated power-spectra plots 
as the compensated plots operate by visually enhancing the high 
wavenumbers at the expense of the low wavenumbers. 



Time (t/to) 


Figure 9. Slopes of the power spectra for the solenoidal (blue 
line) and compressional components (red line) with time for the 
1024® simulation. The inertial range of fc = 20 — 40 as determined 
from the compensated velocity power spectrum was used here to 
plot the power spectral slope at each time interval. The slopes of 
the density spectrum (dashed line) converge to an index value ~ - 
1.2. In comparison, the slopes of the velocity spectrum (solid black 
line) maintain large variations through the entire simulation, but 
have an average value of about -2.0 in the time range 1.0-2.0 in 
agreement with Burgers type turbulence. 


who obtained a slope of -1.95 for the total velocity, -2.02 for 
the compressional components, and -1.92 for the solenoidal 
components. 

The small difference between the solenoidal and com¬ 
pressible spectra suggests that physically, compressional 
shocks may be slightly more efficient at dissipating energy 
compared to solenoidal shocks and hence the slightly steeper 
power-spectrum of the compressible modes. This conclusion 
can be investigated in further detail by examining the in¬ 
terplay between the development of the compressional and 
solenoidal components with time. We do so in Figure[2] where 
we plot the values of the spectral index between fc = 20 — 40 
determined from the power spectra plots for both the den¬ 
sity and velocity spectra with time over the entire simn- 
lation run-time. After some initial spurious values at early 
times before the turbulence saturates, the index of the den¬ 
sity power-spectrum (dotted line) smoothly builds up to 
maintain a relatively stable average value of ~ —1.2 after 
time 1. 0 fp until the end of the simulation. This is in con¬ 
trast to lKritsuk et al.l (l2007l ) who found their density power- 
spectrum exhibited strong variations on short timescales. 
However, in comparison our time evolution of the veloc¬ 
ity power spectrum index is highly variable throughout the 
simulation. The solenoidal (blue) and compressional (red) 
component indexes are seen to sometimes vary by a large 
amount, yet the average spectral index fluctuates around 
a value of -2.0 during the saturated phase. This suggests 
that the nature of our real-space actively driven turbulence 
is highly variable in terms of velocity, yet retains an aver¬ 
age value of -2 representing Burgers strong shock type tur¬ 
bulence. Our spectral slope steepening is due new outflow 


events where we see a strong event can easily change the 
spectral slope. Therefore, when determining power spectra 
of velocity fields with locally forced turbulent simulations, 
one should ensure to average over time otherwise the system 
at an instantaneous time may provide misleading results. 
This is due to the fact that the spectral index of velocity 
power spectra measured in observations could be changed 
a lot by localised strong outflows, supernova explosions, or 
any other energetic point-source driving sources. 

3.5 Ratio of the solenoidal to compressional 
velocity components 

Another method of quantifying the nature of our veloc¬ 
ity field is measuring the relative strengths between the 
solenoid al to compressional co mpo nents, as origina l ly de - 
fined in iFederrath et ^ ll2008l) and iFederrath et al.l (120101 ) 
as 

Vi, = + (1 - ovl = gSi, + (1 - 20 ( 10 ) 

where a Helmholtz decomposition of the projection opera¬ 
tor, Vjj, can represent a purely solenoidal field if ^ = 1, or a 
purely compressible field if = 0, where the projection op¬ 
erators for the fully solenoidal and fully compressive modes 
are Vtj = ^ij — kikj/k^, and = kikj/k^, respectively, 
with Sij representing the Kronecker delta. 

Although the authors defined it with respect to their 
forcing function, we define it based on the resulting turbu¬ 
lent velocity field at each output dump during the simu¬ 
lation. A value of ((=0.0 would represent purely compress¬ 
ible turbulence, and a value of C=1'0 would represent purely 
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Figure 10. Evolution of the relative strengths of the compres- 
sional to solenoidal components, referred to as the value, during 
the 1024^ simulation. We obtain an average of 0.614 between 
a time of 1.0 to 2.0 to indicating that solenoidal modes slightly 
dominate throughout the simulation. It is close to the expected 
2/3 value based on the degree of freedom interpretation. 
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solenoidal turbulence. We plot our ( evolution in Figure fTOl 
Just after time Oto when a few outflows have been intro¬ 
duced to the grid but have not yet interacted, the ( value 
is 0.535 showing each outflow in our model introduces al¬ 
most equal amounts of compressional and solenoidal modes, 
although slightly skewed towards solenoidal modes. After 
the outflows have interacted and the turbulence has become 
fully developed, we measure the average ( to be 0.614 be¬ 
tween time 1.0-2.0 to. This value indicates solenoidal modes 
have increased and slightly dominate in our fully converged 
saturated turbulent velocity field. Interestingly it is close to 
the expected value of 2/3 when we consider the degrees of 
freedom approach where statistically the system will try to 
maintain an equilibrium of 2/3 solenoidal and 1/3 compres¬ 
sional modes. This arises because the compressible mode is 
longitudinal where there is only one degree of freedom (prop¬ 
agation direction). Whereas the solenoidal mode is trans¬ 
verse having two degrees of freedom (two perpendicular 
directions with respect to the propagation direction). The 
slight discrepancy may be due to the fact that a constant 
strong compressional driving is biasing the natural balance. 


4 CONCLUSIONS 

In this paper we study the power spectra of outflow driven 
turbulence through high-resolution isothermal simulations 
using the tvd code. Unlike the alternative method of 
Fourier-space driven turbulence that drives the turbulence 
on global scales, in our model the turbulence is driven lo¬ 
cally in real-space by a se ries of spherical outflows. We use 
the parameter scheme of iMatzneil (l2007l ) in order to ap¬ 
proximate the physical values of a typical star-formation 
cloud such as NGC 1333. The resulting turbulence saturates 
with an average RMS Mach number of ~2.5, and the time- 
averaged power spectrum as well as the relative strengths 
of the solenoidal to compressional components are investi¬ 
gated. 

An index of -1.2 is measured for our time-averaged den¬ 
sity power spectrum. It is a lower value of density power 
spectral slope than found in subsonic compressible flows. 
An explanation is the fact that strong shocks gather the gas 


together into a network of thin dense sheets and filaments. 
Then, in Fourier space, this translates to more power within 
the higher-frequency wavenumbers. Conversely the large vol¬ 
umes of low-density cavities lead to a depletion of lower- 
frequency wavenumbers. Both effects lead to a shallower 
powe r-law slope when viewing in Fourier space (iKim fc Rvul 
l2005li . 

In contrast, our total velocity spectral index is mea¬ 
sured as ^2.0, comparable to Burgers turbulence . This value 
is consistent with the recent work of iFederrathI (l2013f) who 
performed higher Mach number {M — 17) Fourier driven 
turbulence simulations with different total power. This sug¬ 
gests the ability to obtain Burgers shock dominated turbu¬ 
lence through velocity spectra is quite robust despite vari¬ 
ables of Mach number, resolution, and power. 

The nature of power spectra and the Fourier transform 
allow us to readily separate the velocity field into compres¬ 
sional (colliding shocks), and solenoidal (rotational motions) 
components. We see that on average, the slope of the com¬ 
pressional power spectrum (-2.05) is slightly steeper than 
that of the solenoidal spectrum (-1.98). This can be ex¬ 
plained by a combination of strong curved compressible 
shocks generating solenoidal modes, and compressible modes 
decaying more rapidly than solenoidal modes at the small 
dissipation scales leading to a surplus of solenoidal power at 
small scales, and hence a steeper compressible power spec¬ 
tral slope due to the resulting lack of compressible high- 
frequency wavenumbers. The same result can be seen in the 
plot of relative strengths between the solenoidal to com¬ 
pressional components in Figure [10] where we obtain an av¬ 
erage ( value of 0.614, showing solenoidal modes slightly 
dominate after each outflow introduces almost an equal 
amount of solenoidal to compressional mode components 
at early times. The Fourier space driven turbulence models 
of iFederrathI (l2013li also show a steeper compressional than 
solenoidal spectrum. Statistically, compressional modes pos¬ 
sess 1 degree of freedom, whereas solenoidal modes possess 2 
degrees of freedom. The system tries to maintain equilibrium 
of 1/3 compressible and 2/3 solenoidal. This is consistent 
with our results. 

The average density weighted velocity power spectrum 
of our model is measured with a spectral index of 
-1.6. This value is slightly smaller than the Kolmogorov -5/3 
scaling found in previous studies at similar resolution. If un¬ 
certainties of resolution and the short inertial range fitting 
could be excluded, it may be due to the nature of our real- 
space driving mechanism and the low value of the spectral 
index of the density field. This may suggest power spec¬ 
tra reveal a fundamental difference between Fourier based 
driven turbulence and real-space driven turbulence. Fourier 
based driving could be considered as large scale driving; it 
stirs the entire computational domain equally. Then the tur¬ 
bulence must decay transferring energy from large scale to 
small scale. In contrast real-space driving only stirs the com¬ 
putational domain in localised regions. Energy needs to be 
transferred up from small scale to large scale, against the 
natural direction of turbulence decay resulting in more en¬ 
ergy remaining locked at small scales. As the spherical out¬ 
flows in our model expand, they transfer power to a larger 
scale than their initial injection scale. This transfer of power 
from a small scale to a larger scale is not due to turbulence, 
but the nature of our real-space driving model. 
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We tentatively propose that the density weighted ve¬ 
locity power spectrum is not a universal power spec¬ 

trum. Once Mach number reaches a certain value, the index 
of total velocity is -2.0, even if driven with much greater 
Mach number. In contrast, the index of density power is vari¬ 
able. So a combination of density and velocity power would 
still be variable with Mach number. We are currently inves¬ 
tigating how the spectrum of depends on the Mach 

numbers for a future publication. 

Returning to the consideration of NGC 1333, our simu¬ 
lation results suggest that in an outflow-driven environment 
the solenoidal velocity component may be stronger than the 
compressible velocity component, and the spectral slope of 
density power spectrum may be closely related to the global 
Mach number of the medium. 
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